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Abstract 

In Vuik and Ryan (2014) we studied the use of troubled-cell indicators for dis¬ 
continuity detection in nonlinear hyperbolic partial differential equations and intro¬ 
duced a new multiwavelet technique to detect troubled cells. We found that these 
methods perform well as long as a suitable, problem-dependent parameter is cho¬ 
sen. This parameter is used in a threshold which decides whether or not to detect 
an element as a troubled cell. Until now, these parameters could not be chosen 
automatically. The choice of the parameter has impact on the approximation: it de¬ 
termines the strictness of the troubled-cell indicator. An inappropriate choice of the 
parameter will result in detection (and limiting) of too few or too many elements. 

The optimal parameter is chosen such that the minimal number of troubled cells is 
detected and the resulting approximation is free of spurious oscillations. 

In this paper we will see that for each troubled-cell indicator the sudden in¬ 
crease or decrease of the indicator value with respect to the neighboring values is 
important for detection. Indication basically reduces to detecting the outliers of a 
vector (one dimension) or matrix (two dimensions). This is done using Tukey’s box- 
plot approach to detect which coefficients in a vector are straying far beyond others 
(Tukey, 1977). We provide an algorithm that can be applied to various troubled-cell 
indication variables. Using this technique the problem-dependent parameter that 
the original indicator requires is no longer necessary as the parameter will be chosen 
automatically. 
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1 Introduction 

In [Snj , we studied the use of troubled-cell indicators for discontinuity detection in non- 
linear hyperbolic partial differential equations and introduced a new multiwavelet tech¬ 
nique to detect troubled cells. We compared the troubled-cell indicator of Qiu et al. 
using Harten’s subcell resolution [m [28] , the KXRCF shock detector [23| and the mul¬ 
tiwavelet troubled-cell indicator [39]. We found that these methods perform well as long 
as a suitable, problem-dependent parameter is chosen, which was also observed in [28]. 
This parameter is used in a threshold which decides whether or not to detect an element 
as a troubled cell. Until now, these parameters could not be chosen automatically such 
that the indicator works well in a variety of situations m- Similarly, a parameter is 
required for adaptive mesh refinement HU- Here, the used threshold parameter does 
depend on the discretization and the order of accuracy. This parameter is always chosen 
in the same way, similar to the KXRCF method [23]. 

The choice of the parameter has impact on the approximation: it determines the 
strictness of the troubled-cell indicator. An inappropriate choice of the parameter will 
result in detection (and limiting) of too few or too many elements. Detection of too 
few elements leads to spurious oscillations, since not enough elements are limited. If 
too many elements are detected, then the limiter is applied too often, and therefore 
the method is more costly and the approximation smooths out after a long time. The 
optimal parameter is chosen such that the minimal number of troubled cells is detected 
and the resulting approximation is free of non-physical spurious oscillations. In general, 
many tests are required to obtain this optimal parameter for each problem [25] 02] • 

In this paper we will see that for each troubled-cell indicator the sudden increase 
or decrease of the indicator value with respect to the neighboring values is important 
for detection. Indication basically reduces to detecting the outliers of a vector (one 
dimension) or matrix (two dimensions). This is done using Tukey’s boxplot approach 
to detect which coefficients in a vector are straying far beyond others m, which is 
commonly used in statistical analysis uniEH]. This method is designed in such a way 
that only a few ’false positives’ are found if the data are well behaved (i.e., Gaussian 
USD- Another advantage of this method is that it is not necessary to specify the number 
of possible outliers in advance. This is in contrast to many standard outlier-detection 
techniques which require a statement of the exact or the maximum number of outliers 
that may be present m- 

We provide an algorithm that can be applied to various troubled-cell indication vari¬ 
ables. Using this technique the problem-dependent parameter that the original indicator 
requires is no longer necessary as the parameter will be chosen automatically. 

The numerical results in this paper are computed using the discontinuous Galerkin 
(DG) method [71[6l[5l[8] together with a third-order strong stability-preserving Runge- 
Kutta time-stepping scheme m- We apply either the original troubled-cell indicators 
(with an optimal parameter), or the outlier-detection technique in combination with the 
indication variable. In that way, the performance of the new technique can be easily com¬ 
pared to the current state-of-the-art methods. We will apply the techniques to various 
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test problems in one and two dimensions. Here, we investigate the modified multiwavelet 
troubled-cell indicator m, the KXRCF shock detector [23] and the minmod-based TVB 
indicator [7] in more detail. The moment limiter is used in the detected troubled cells 
[22], but other limiting techniques can be used. 

The outline of this paper is as follows: in ^ we present the relevant background 
information on discontinuous Galerkin methods, troubled-cell indicators and the moment 
limiter. In ^we introduce our new outlier-detection algorithm. The effectivity of this 
new method compared with the corresponding parameter-using troubled-cell indicators 
is presented in 0 for standard numerical examples. The computational costs of our 
algorithm are addressed in ^ We conclude with a discussion of our method and future 
work in ^ 

2 Background 

This section contains some background information about the discontinuous Galerkin 
(DG) method [7l[6l[5l[8|, as well as the theory behind troubled-cell indicators [SaClEHl 
Esmo] and the moment limiter [22] • This information will be used to apply our new 
outlier-detection scheme in numerical examples. 

2.1 Discontinuous Galerkin method 

We briefly explain the DG method using the initial-value problem 

utAf{u)x = 0, X G [-1,1], t > 0, , . 

u{x, 0) = u^{x), X G [—1,1], 

where u = u{x,t), and f{u) describes the flux function. 

Discretization in space is obtained by dividing [—1,1] into 2^ elements (used in the 
multiwavelet expansion, ^2.2.11) . defined as 

Ij = [Xj_i,x.^), j = 0,... ,2^ - 1. 

The choice for half-open intervals follows from the paper of Archibald et al. [3]. Different 
choices are available in the literature, for example closed intervals (Hovhannisyan et al. 
m), or open intervals (Gerhard et al. [11]). 

The approximation space that we use on each element is Vh{Ij) = {x G P^(/j)}, where 
is the space of polynomials of degree k. In order to take advantage of the multiwavelet 
properties, the basis for P^ is constructed using the scaled Legendre polynomials, which 
are defined as 

where is the Legendre polynomial of degree I, I — t),... ,k. 


( 2 ) 
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The weak form of the PDE in problem ([T]) is obtained by multiplying the equation by 
a test function v G Vh{Ij) and integrating over element Ij. Using integration by parts, 
this yields 


/ Utvdx = / 

Jp Jp 


f{u)v^dx + 1 - A-lit’ 


2 J-o 




(3) 


where fj±i denote the flux values through boundaries Xj±ii 2 - These are approximated 
using the local Lax Friedrichs flux [2S]: 



1 

2 




where we use that / is convex, such that 

a-i = max(|/'(u" i)|, |/'(ut i)|). 

•^2 J 2 J 2 

The third-order strong stability-preserving Runge-Kutta scheme m is used for time 
evolution. Note that this is only a choice, and that other time-stepping schemes are also 
possible [niEiiissj. 


2.2 Troubled-cell indicators 

In this section, various troubled-cell indicators are described on which our new outlier- 
detection algorithm will be tested. In particular, we will investigate the modified mul¬ 
tiwavelet troubled-cell indicator 133110], the KXRCF indicator [23], and the minmod- 
based TVB indicator [Tj. In earlier papers, Harten’s subcell resolution idea mi was 
used for indication [23 EH]. However, this method was unstable in several numerical 
experiments [33] ? therefore, we will not investigate this method here. 

The outlier-detection algorithm will require that we pass in a vector of troubled-cell 
indication variables and therefore we provide the vector form as well. 


2.2.1 Modified multiwavelet troubled-cell indicator 

In [39], a multiwavelet troubled-cell indicator was constructed, which was modified in m 
by using multiwavelet coefficients instead of multiwavelet contributions. In this section, 
we repeat the important definitions. Here, we only investigate the domains [—1,1] (one 
dimension) and [—1,1] x [—1,1] (two dimensions). The corresponding definitions can be 
easily extended to general domains in one and two dimensions [39] . 

A global DG approximation of degree k on 2'^ elements in [—1,1] can be written as 

2^-1 k 

Uh{x) = 2"t X e [-1,1], 

j=0 i=0 

where the scaling functions (jf^- are defined as 

4>l{x) = 2^l^4>t{2^{x + 1) - 2j - 1), ^ = 0,..., fc, j = 0,... , 2^^ - 1, 
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and are the scaled Legendre polynomials as in equation ([2]). 

The corresponding multiwavelet decomposition of the DG approximation can be 
written as 

k n-l 2^-1 k 

Uh{x) = '^s%(l)e{x) + Y^ 

i=0 ni=0 j=0 i=0 

where are the scaling-function coefficients belonging to u^, the multiwavelets on 
higher levels are defined as + 1) — 2j — 1), and djj are the cor¬ 

responding multiwavelet coefficients 13 EH], which are determined using the orthogonal 
projection of the DG approximation onto the multiwavelet basis: 

^ij — '0^)L2([-l+2-"^+lj,-l+2-"^+bJ+l)])• 

In practice, these coefficients are efficiently computed using the quadrature mirror filter 
coefficients 13 Ej. The multiwavelets have been developed by Alpert [I], and are also 
explained in m- 

As we have seen in m, the coefficients on level n — 1 are strongly related to the 
inter-element jumps in (the derivatives of) the DG approximation. The multiwavelet 
coefficients on level n — 1 of the decomposition equal 


<•7'=E c., ■ («r 

(^^■+ 1 / 2 ) ^(^2j+l/2)) ’ 

(4a) 

m=0 



o(—n+l)m 

m! 

f x'^^i{x)dx, 

Jo 

(4b) 


where £ = 0,..., /c, j = 0,..., 2^“^ — 1, and is the mth derivative of Note 
that this only includes half of the element-boundary jumps. In order to also compute 
the rest of the jumps, we use the same renumbering technique as was proposed in gnj. 
This gives rise to 2'^ coefficients for level n — 1, in the following denoted by d^- ^, £ = 
0,...,A:,j = 0,...,2--l. 

In general, the DG approximation is discontinuous at element boundaries. Therefore, 
the multiwavelet coefficients are usually not exactly equal to zero. However, when the 
solution is sufficiently smooth, then the element-boundary jumps in the approximation 
and its derivatives will be noticeably smaller than when a discontinuity in (one of the 
derivatives of) the solution is present due to the cancellation property of multiwavelets 
gHj. The multiwavelet coefficients ^ are used to detect troubled cells when 

^ = 0,... ,2^ - 1}, G G [0,1]. (5) 

Since coefficient contains information about the jump in (derivatives of) the DG 
approximation at Xj+ 1 / 2 , elements Ij and Ij^i are limited if d^J^ satisfies inequality ([5]). 
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Note that only the multiwavelet coefficients on level n — 1 are used for indication. 
Therefore, a twoscale representation of the DG approximation would suffice, and the 
number of elements in the domain might be even (instead of the restriction to a power 
of two). However, since this would change the definitions, we have chosen to use 2^ 
elements in this work. 


The boundary of which is maximal, is assumed to be the location where the 

strongest shock occurs. If C = 1 , then no element will be detected, and the smaller C is, 
the more elements will be limited. In this way, the value of (7 is a useful tool to prescribe 
the strictness of the limiter. In general, it is hard to choose a sufficient value for C. For 
each problem, several tests should be done in order to obtain an optimal parameter [39] . 


In order to remove the problem-dependent parameter C that occurs in this indicator, 
we propose to use the multiwavelet coefficients in our outlier-detection algorithm (Q. 
The indication vector is defined as D = (<^^ 0^5 • • • 5 <^^ 2 ^- 1 )^- 

In two dimensions, the relations for the multiwavelet coefficients on level n — 1 follow 
naturally from the one-dimensional coefficients: 


7 a,n-l _ 


k 

E C,.,/ " 


^ 2 -^ y 

/ -J I 

mx=0 y^.n-k 


ifd^yuh^ + ^ d^yuh^ _ ^ \ 


dyy‘ 




k 


( d^^Uh ^ , y) - 7""^^ (a;- i,y) ) .\y)dy, 


E E 

ma;=0 my=0 




.(■ 


Qrrix gruy 


'^^x \ ddy'^^y 


where £ = ( 4 , 4 )^, j = and n = i = 0, ..., 2 ""^“^ - 1 and j = 

0,..., 2'^y~^ — 1, and ^ is defined as in equation f|4b|) . Using the renumbering technique 
as mentioned before m, we include the second half of the boundaries in the multiwavelet 
direction, and find coefficients (j = 0 ,..., 2^y — 1), {i — 0,, 2^^ — 1) 

^ (i = 0,...,2-^-1, j = 0,...,2-^-1). 

Note that these relations indeed confirm the observations that the ol mode detects 
discontinuities in the y-, the /3 mode in the x-, and the 7 mode in the x 7 /-direction, as 
was stated in [2H] and seen in [59] . 

In two dimensions, the approach of inequality ([5]) is applied for each mode separately. 
In the a mode, we take the coefficients with index £ — (0, k)^ for indication. In the (3 
mode, the indices £ = (fc, 0 )^ are used, and for 7 we take £ = (fc, k)^. 


and ^ 


Using the directions of each mode, the one-dimensional outlier-detection algorithm is 
applied to the a-mode vectors for each x, and to the /3-mode vectors for each y. We have 
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found that detection on the 7 mode selects too many elements. Therefore, this mode 
is not used in the outlier-detection scheme. We apply outlier detection to the following 
vectors: 


a : D, 


( ra,n—l Ta,n—1 \ • _ p. 


n _ f 7/3,n—1 7/3,n—1 




, J = 0,...,2 


2nx-l _ 

Tiy 1 _ ^ 


We note that the multiwavelet indicator is equivalent to comparing the DG approxima¬ 
tion over two different spatial meshes. An alternative strategy would be to compare the 
approximation at two different levels in time as in Dumbser et al. [9]. 


2.2.2 KXRCF indicator 


The shock-detection technique by Krivodonova et al. [23] uses inffow boundaries to 
detect troubled cells. The detector considers the jump in across the inffow edges of 
Ij and examines 


7 = 



Here, dij is the inffow boundary and Uh\i^. is the DG approximation in the neighbor 
of Ij on the side of d/ 7 . The indicator is normalized to 

J 


7- 


Ioi-iuh\ij -Uh\i )ds 


h^\dl-\\\u,\. 


J =0,...,2”-l. 


( 6 ) 


Here, h is the radius of the circumscribed circle in Ij, and the norm is based on the av¬ 
erage in one dimension and the maximum norm in quadrature points in two dimensions. 

Near a discontinuity Ij ^ 00 , whereas Xj ^ 0 if /i ^ 0 or fc ^ 00 in smooth-solution 
regions. In [23], the threshold value is taken equal to 1, such that element Ij is detected 
as troubled if Ij > 1, and in that case the limiter is applied in Ij. Note that this 
threshold parameter is chosen arbitrarily: the value 1 does not necessarily follow from 
the theory. 


In order to remove this parameter, the outlier-detection mechanism was tested on 
the vector D = (Iq, ... ,l 2 n-i)^. However, it turned out that the original discontinuity 
detector without normalization is more suitable, such that the jump across the interfaces 
is used in the indicator: we take D = (Xq, ... ,l 2 ^-i)^ for detection. 

In two dimensions, a matrix D = {Ijj} is found. Here, the one-dimensional outlier- 
detection approach is applied in the x- and ^/-direction separately (row and column 
wise). 
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2.2.3 Minmod-based TVB indicator 


In this section, the minmod-based TVB indicator of Cockburn and Shu will be explained 
id- For each element Ij,j = 0,... , 2^ — 1, the element-boundary approximations are 
split into 


u 1 = Uj -h Uj and = Uj — Uj, 

^'2 


(7) 


where 

k k 

£=1 i=l 

Element Ij is detected as troubled if either Uj or Uj is modified by the functions 

-(mod) _ _ Uj,Uj — Uj-i), = m{Uj,Ujj^i — Uj,Uj — Uj-i), (8) 

where the TVB-modified minmod function is defined as 

^ . \ _ { Gil’) if |«-i| < MAx^, 

m(ai,..., Uq) , Uq), otherwise, 

in contrast with the standard minmod function 

, _ r s • \aj\, if sign(ai) = • • • = sign(ag) = s, 

1 0 , otherwise. 


(9) 


Note that the parameter M is difficult to tune, and hardly any difference is found when 
M ranges from 1 to 100, m- We use the minmod-based TVB indicator for detection 
and then apply a chosen limiter in the detected troubled cells. 

For systems of equations, characteristic field decompositions are required [H]. The 
corresponding eigenvector matrix is computed using Roe averages 0130]. 


Instead of using the parameter M, we will apply the outlier-detection algorithm. DG 
(1) (k) 

coefficients Uj \ ... ,Uj usually differ substantially from their neighbors when Ij belongs 
to a discontinuous region. This means that we use the vectors Di = {uq, ... ,U 2 ^-i)^ 
and D 2 = (f/Q? • • • ^U 2 n-i)^ in the outlier-detection technique and detect element Ij as 
troubled if either Uj or Uj is detected as an outlier. 

For two-dimensional systems, the procedure for has been explained in [8]. The 
indicator uses solution derivatives (e.g. DG coefficients) for detection. We use Q^, 
which means that more ’cross-product’ coefficients exist (for example, for k = 1: 
However, using Biswas’s reasoning [3], we do not use these coefficients for detection, 
since they have a lesser effect on the numerical approximation than either or • 
The minmod-based TVB indicator in two dimensions resembles the two-dimensional 
moment limiter [22]. The difference between the two approaches is that the moment 
limiter uses forward and backward differences of lower derivatives, whereas the minmod- 
based indicator uses a finite-difference approach on the element averages. 
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In our numerical examples, we will focus on the case A: = 1, and use the DG coeffi¬ 
cients and for detection. Outlier detection will be applied to the vectors 

— \^0J ’ • • •: 


y 

II 

o 

- 1, 

in the x-direction and 







y 

, i = 0, 

2n. _ 


in the ^-direction. In this way, it is possible to detect discontinuities in different directions 
as we will see in ^4.2[ 


2.3 Moment limiter 

In the detected troubled cells, a limiter is applied. The limiting technique that we use 
in this paper is the moment limiter [22]. This is only a choice - other limiters are also 
possible. 

The moment limiter reduces the DG approximation to a low order in discontinuous 
regions, and maintains a high order if the approximation is smooth enough. Although 
the limiter has its own mechanism to control which regions should be limited, we will 
apply troubled-cell indicators as a switch to control where the limiter is applied. This is 
to prevent limiting smooth extrema. 

The moment limiter limits DG coefficients, starting at the highest level k. For each 
element Ij,j = 0 ,..., 2^ — 1, the limited value of coefficient equals 



with Pk — {\/k — l/2)/(^/c + l/2) and using the minmod function (equation ([9])). If 
5^-^^ = then the limiting procedure is cut off for this element Ij. If not, then 
is limited using the same procedure, continuing until is limited, or stopping the first 
time for some ^ = A; — 1,..., 1. 

(i) 

For systems of equations the limiter is applied to the characteristic variables w^- ^ = 

—1 (t) 

R Aij. Due to this approach it is possible that negative values for density, pressure 
or energy are found. In that case, all higher-order coefficients are set equal to zero, and 
is limited using equation m- If negative values are still found, then the linear 
coefficient is also set equal to zero. 

In two dimensions, the moment limiter uses the neighboring elements both in the x-, 
and in the ^-direction [22]. 


3 Troubled-cell indication using outlier detection 

In this section, an outlier-detection algorithm is proposed to detect outliers in a vector. 
This technique will be applied to the troubled-cell indicators given in ^2.21 
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In order to detect outliers we use a boxplot mechanism that is often applied in 
statistics uni EH!, and described by Tukey 137j. Important properties of this method are 
that only a few ’false positives’ are found if the data are well behaved (i.e., Gaussian 
USD, and that it is not necessary to specify the number of possible outliers in advance. 
This is in contrast to many standard outlier-detection techniques which require to state 
the exact or the maximum number of outliers that may be present jl5j . 

Here, we use a vector d = {do, ..., , of which outliers (suddenly changing coef¬ 

ficients with respect to neighbors) should be detected. A general outline of the outlier- 
detection algorithm that we use is provided below. In the following we discuss the 
details. 

Algorithm 1 Outlier-detection algorithm. 

Send in a suitable troubled-cell indication vector d. 

Sort d to obtain d^. 

Compute the quartiles of d^. 

Construct the outer fences. 

Determine the outliers. 


3.1 Quartiles 

Quartiles separate the data into four equal groups [29] . The values of Qi, Q 2 (the 
median), and Qs provide useful information about the structure of d. As a preparation, 
it is convenient to sort d, such that we obtain the vector d^: 

d^ = {do, dl,..., dp^)^, where, do — — • • • — 

The median of d is defined as the numerical value that separates the higher half of the 
vector from the lower half m- It equals 

5 if ^ is even, 

1 {dlN-l)/2 + ^y+i)/2) , if ^ is odd. 

The median is also called the second quartile of the vector d. 

The first quartile is defined as the value below which 25% of the data fall, and is 
denoted by Qi. Similarly, the third quartile, Q 3 , equals the value that splits off the 
lowest 75% of the data from the highest 25% [29|. Many different definitions of the first 
and third quartiles are used. In this work we apply Tukey’s definition (definition 6 in 

uni): 

Qi ^ - 9)dj-i+gdj, ( 11 ) 

where [{N -h 4)/2]/2 — j + and [x] denotes the largest integer that does not exceed 
X. Note that ^ = 0 or ^ = 1/2. The third quartile Qs is then computed symmetrically 
using the upper end of the vector d^. 

In practice, we will always use a vector with N + 1 = 4r coefficients, where r G N. 
In that case, Qi = (d^_i + d^)/2 and Qs = (<^ 37—1 + <^ 3 r)/ 2 - 
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3.2 Fences and outlier detection 

We have already seen that the values of the quartiles provide useful information about 
the structure of d^. However, this is not enough to define outliers in the vector. Outliers 
are the coefficients in the vector that are straying far out beyond the others. In order 
to pick out certain coefficients as outliers, inner and outer fences are constructed, which 
were originally defined by Tukey m- The inner fences are equal to [Qi — 1.5{Qs — 
Qi)^Q3 + 1-5((53 — Qi)] (coefficients outside this interval are called soft outliers). When 
the data are normally distributed, only 0.7% of the data set is seen as a soft outlier 
(asymptotically) [15]. The value 1.5 is referred to as the whisker length of the boxplot. 

The outer fences of a vector are [Qi — 3{Qs — Qi),Q3 + — Qi)] (coefficients 

outside are called extreme outliers). The coverage for this whisker length is 99.9998%, 
such that only 0.0002% of the data in a normally distributed vector is detected as an 
extreme outlier (asymptotically) [15]. The choices of the whisker lengths (1.5 and 3) were 
proposed by Tukey and are commonly used in the literature [IQI[I51[I91[201[3I1[32]. 
We will use the extreme outliers to detect troubled cells, since then very outstanding 
coefficients in the vector are selected. Because the data were sorted, the outer fences 
and outliers can easily be determined. 

3.3 Application to troubled-cell indication variables 

In this section, we will explain the application of outlier detection to troubled-cell indica¬ 
tion variables in one dimension. The corresponding indication vectors for each troubled¬ 
cell indicator were given in ^2.2[ In this section, we have seen that all described indicators 
attach a value to each element of the domain (multiwavelet coefficient, jump across inflow 
boundary, or approximation at boundaries). Discontinuous regions usually correspond 
to the locations where the indicator value suddenly increases or decreases with respect 
to the neighboring values. This means that indication basically reduces to detecting the 
outliers of a vector with troubled-cell indication values. By applying the new outlier- 
detection technique, the threshold to be an extreme outlier is fixed, and the indicator 
no longer depends on problem-dependent parameters. 

When an approximation contains several discontinuous regions, outlier detection ap¬ 
plied to the global vector D will only select the strongest discontinuities. In order to 
also take into account the weaker discontinuities and the local structure of the approxi¬ 
mation, the vector D will be split into local vectors of fixed length. For each subvector 
the outlier-detection mechanism is applied. In the local approach we ignore the detected 
coefficients in the left half of the local region if they are not detected with respect to 
the left-neighboring vector, and similarly the detected coefficients in the right half of the 
local region are tested. In this way the spatial information can still be used. 

The outlier-detection algorithm executes the steps as provided in Algorithmic Below 
we explain this in more detail. 

Since the global vector D consists of 2'^ coefficients, we propose to split D into 
local vectors of length 2^, where p G {2,... ,n}. Each local vector is then sorted. For 
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Algorithm 2 Outlier-detection algorithm using local vectors. 

Send in a suitable troubled-cell indication vector D. 

Split this vector into local vectors, d. 
for all local vectors do 
Sort d to obtain d^. 

Compute Qi and Qs using definition ffTTl) . 

Detect dj in the smallest 25% of d^ if dj < Qi — 3{Qs — Qi), and dj in the biggest 
25% of d^ if > Qs + 3(Q3 - Qi). 

end for 

Ignore the detected outliers in the left half of the local region when they are not 
detected with respect to the left-neighboring vector, and similarly test the detected 
coefficients in the right half of the local region. 


convenience, we denote the sorted local vector of coefficients by d^ = (dg, df,... d^), 
where N — 2^ — 1. By definition this vector has the following 25th and 75th percentiles 
(see equation (ED): 


Qi = 


— 1 % 


Qs = 


^3-2p-2_x ' ^3.2p-2 


Next we compute outer fences. Outliers are determined by comparing the smallest vector 
values with Qi — 3{Qs — Qi) and the biggest components with — Qi). For 

the smallest values we start with testing whether dg < Qi — 3 (Q 3 — Qi). If dg is not an 
outlier, then there are no other outliers, since dj > d^ > Qi — 3{Qs — Qi), j = 0 ,..., A. 
If dg is an outlier, then we test df, etcetera. By construction, Qi — 3{Qs — Qi) < Qi, 
such that the only possibilities for low outliers are dg,..., d^p_ 2_2 ( 2 ^“^ — 1 coefficients). 
This means that at most dg, df,..., d 2 p- 2_2 should be tested. 

Similarly we test df^ and (possibly) ..., d^ 2 ^- 2 ^]^ against (^3 -h 3{Q^ — Qi) 

(depending on the outcome). Also here, at most 2 ^“^ — 1 coefficients should be tested. 

Finally, the detected outliers in the left half of the local vector are compared with 
the fences of the left-neighboring vector, and the outliers in the right half are compared 
with the right-neighboring fences. 

Considering the number of elements in each local vector, it should be noticed that 
p = 3 results in 8 coefficients per vector, which is too few to find a boxplot which 
is meaningful. Using p = 4 (16 coefficients per vector) means that at maximum six 
outliers can be detected per local vector. Therefore, the maximum number of possible 
outliers in D equals 2^“^ •6 = 3* 2'^~^. If we take more coefficients per local vector, 
for example p = 5 (32 coefficients per vector), then the ’stencil’ is too big to extract all 
local information of the approximation. Therefore, we propose to use 16 coefficients per 
local vector (p = 4), which worked well in all test cases we performed. 

In two dimensions, the one-dimensional algorithm is applied in the x- and p-direction 
separately. The corresponding troubled-cell indication vectors were given in ^ 2.21 
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4 Numerical results 

In this section, the original troubled-cell indicators are compared with the new outlier- 
detection approaches. This is done for the modified multiwavelet troubled-cell indicator 
of Vuik and Ryan m, the KXRCF indicator [23] , and the minmod-based TVB indicator 
[Tj. We computed the results using k = 1,2,3. In this paper, we only present the case 
k = 2. 

The results for the one-dimensional test cases are presented using time-history plots 
of detected troubled cells, which is commonly done [391 HOl H3] . 

4.1 One-dimensional tests 

The test cases in one dimension include one continuous example using the Euler equations 
on [—1,1] with initial conditions po{x) = 1 -h 0.5 sin(107rx), uo{x) = 1, po{x) = 1, and 
periodic boundary conditions. The solution at final time T = 2 is given by p(x, 2) = 
Po{x). Using this example, we can validate our algorithm: since no discontinuities are 
present, no element should be detected. Indeed, the original troubled-cell indicators do 
detect certain elements (chosen parameters are reasonable, and commonly used [71 ESI 
[SS]). This is depicted in Figure [U in which the detected troubled cells using the original 
indicators are visualized. These so-called time-history plots show which elements are 
detected in space for each time step. 

The application of the outlier-detection algorithm together with the troubled-cell 
indication vectors does not select any element, which is the desirable result. 



-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 


(a) Multiwavelets, C = 0.5 (b) KXRCF, threshold 1 (c) Minmod, M = 10 

Figure 1: Detected troubled cells for p(x,2) = 1 -h 0.5 sin(107rx), 128 elements, fc = 2, 
using original troubled-cell indicators. Corresponding outlier-detection approaches do 
not detect any element. 

The standard numerical examples for the Euler equations are also investigated. 

Below the results and comparisons are shown using four different sets of initial con¬ 
ditions: the shock tubes of Sod [35] (Figure E|) and Lax [23] (Figure [3D, the blast-wave 
problem [H] (Figure ED, and the Shu-Osher problem [33] (Figure (ID . We omit the details 
of these test problems and refer to |39| for more information on initial conditions and 
boundary conditions. We apply the indication technique to density for the modified 
multiwavelet indicator, density and energy for KXRCF, and the characteristic variables 
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for the minmod-based TVB indicator, as has been done by Qiu et al. [2H]- The first row 
of each figure consists of time-history plots of detected troubled cells using the original 
indicators. The second row belongs to the outlier-detected troubled cells. The corre¬ 
sponding approximations at the final times are given in the third and fourth row. Note 
that these results are computed using the moment limiter in the detected troubled cells. 
A different choice for the limiter will result in different approximations. In all figures we 
take k — 2, and similar results were found for A: = 1, 3. 

Note that the original troubled-cell indicators are applied using the optimal problem- 
dependent parameters as found in [281 [39]. We stress that the outlier-detected results 
are computed without problem-dependent parameters, but with a fixed whisker length 
equal to 3, and with local indication vectors of size 16. 

It turns out that the new outlier-detection approach detects the troubled regions 
very accurately and generally better than the original parameter-using methods for the 
blast-wave and Shu-Osher problem. For the shock tube problems of Sod and Lax, most 
discontinuous regions are selected. Note that the outlier-detection indicators sometimes 
detect jumps in derivatives, as can be seen at the end points of the rarefaction waves. 
The original indicators, however, do not detect these structures. This difference can 
be explained by recalling that the original indicators focus on the actual value of the 
indication variable, whereas the outlier-detection techniques investigate the relative value 
with respect to the neighboring region. A discontinuity in the derivative usually causes 
sudden differences, and therefore these regions are detected as troubled. By applying a 
limiter at these locations, the discontinuity in the derivative is smeared a bit, such that 
at some time steps these elements are not detected. Note that all approximations are 
very accurate and close to the exact solution. 

The most important improvements are found for the blast-wave and Shu-Osher prob¬ 
lem. For the blast waves, the original KXRCF detector and minmod-based TVB indica¬ 
tor detect many elements. However, the new outlier-detection approach combined with 
these detection variables only selects a few of them, thereby still producing very accurate 
results. 

In the Shu-Osher problem (Figure [5]) an initial discontinuity is moving to the right, 
thereby evolving (highly oscillatory) continuous regions and developing new shocks in 
the left side of the domain. 

The first row of the figure consists of time-history plots of detected troubled cells 
using the original indicators. Note that both the multiwavelet indicator with C — 0.01 
and the minmod-based TVB indicator with M — 100 detect the highly-oscillatory region 
as being discontinuous. In this case, the KXRCF indicator gives more accurate results. 
For k — 1 however, the KXRCF indicator only detects the largest discontinuity, and 
neglects the other three shocks in the left side of the plot, which leads to some spurious 
oscillations in the approximation. 

In the second row of the figure, the time-history plots are shown when the indication 
vectors are used in the outlier-detection algorithm. All three indication techniques detect 
the correct regions, and the approximations are as expected (row 3-4 of Figure [5]). Note 
that the results are very close to the exact solution: the outlier-detection algorithm is 
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indeed able to replace the problem-dependent parameters in the original indicators. 

For fc = 1 and fc = 3, the same behavior is found: the new outlier-detection approach 
perfectly selects the discontinuous regions in the domain. 



(a) Original, C = 0.1 



(d) Outlier, multiwavelets 



(g) Original, C = 0.1 




(b) Original, KXRCF 



(e) Outlier, KXRCF value 



(h) Original, KXRCF 



1.8 

1.6 

1.4 

1.2 

1 
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-5 -4 -3 -2 -1 0 1 2 3 4 5 


(c) Original, M = 10 



(f) Outlier, minmod-based TVB 



(i) Original, M = 10 



(j) Outlier, multiwavelets (k) Outlier, KXRCF value (1) Outlier, minmod-based TVB 

Figure 2: Detected troubled cells (row 1 and 2) and approximation at final time T = 2 
(row 3 and 4), shock tube of Sod, k — 2, 128 elements. 
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(a) Original, C = 0.1 (b) Original, KXRCF (c) Original, M = 10 





(d) Outlier, multiwavelets 



(e) Outlier, KXRCF value (f) Outlier, minmod-based TVB 



(g) Original, C = 0.1 



(h) Original, KXRCF (i) Original, M — 10 



(j) Outlier, multiwavelets (k) Outlier, KXRCF value (1) Outlier, minmod-based TVB 

Figure 3: Detected troubled cells (row 1 and 2) and approximation at final time T = 1.3 
(row 3 and 4), shock tube of Lax, k — 2, 128 elements. 
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(a) Original, C = 0.05 (b) Original, KXRCF (c) Original, M = 100 



(d) Outlier, multiwavelets (e) Outlier, KXRCF value (f) Outlier, minmod-based TVB 


1-Exact 

' ] - - - Multiwavelets 



6 

5 

4 

3 

2 

1-Exact 

1---KXRCF 



6 

5 

4 

3 

2 

1 

1 - Exact 

]-Mlnmod 



f 

\l 

f 

\l 

f 

\l 

0 0.2 0.4 0.6 

(g) Original, C - 

0.8 1 

= 0.05 


) 0.2 0.4 0.6 0.8 1 

(h) Original, KXRCF 

°( 

) 0.2 0.4 0.6 

(i) Original, M 

0.8 1 

= 100 

\ -Exact 

- - Multiwavelets 

f 

\l 


6 

5 

3 

2 

0 

1-Exact 

1-KXRCF 

f 

\| 


6 

5 

4 

3 

2 

1 

0 

J - Exact 

]-Mlnmod 

f 

\1 



(j) Outlier, multiwavelets (k) Outlier, KXRCF value (1) Outlier, minmod-based TVB 

Figure 4: Detected troubled cells (row 1 and 2) and approximation at final time T = 
0.038 (row 3 and 4), blast-wave problem, k — 2, 512 elements. 
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(a) Original, C = 0.01 


(b) Original, KXRCF 


(c) Original, M = 100 
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(d) Outlier, multiwavelets 



(g) Original, C = 0.01 



(j) Outlier, multiwavelets 


(e) Outlier, KXRCF value 
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(i) Original, M = 100 
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Figure 5: Detected troubled cells (row 1 and 2) and approximation at final time T = 1.8 
(row 3 and 4), Shu-Osher problem, k = 2, 512 elements. 
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4.2 Two-dimensional test 

In two dimensions, we investigate the double Mach reflection of a strong shock which 
satisfies the two-dimensional Euler equations. Again, the original troubled-cell indicators 
(with optimized parameter) are compared to their outlier-detection approaches. The 
results for fc = 2 can be compared in Figure [6] for the modified multiwavelet troubled¬ 
cell indicator, Figure 0 for the KXRCF shock detector, and Figure [H] for the minmod- 
based TVB indicator (A: = 1 only). The spatial domain is split into 2^ x 2^ rectangular 
elements, such that Ax = Ay = 1/128. In each figure, the left plots are computed 
using the original troubled-cell indicators, and the right plots correspond to the outlier- 
detection approaches. 

As mentioned earlier, the multiwavelet technique is able to distinguish between x- 
and ^/-directed discontinuous regions. This is also the case when outlier detection is 
used. We point out that a sharp detection of the discontinuous region is found. Only a 
few elements outside the discontinuous region are added, which apparently correspond 
to discontinuities in derivatives (since jumps in multiwavelet coefficients are found). The 
approximations at the final time are comparable to the results using the original modified 
multiwavelet troubled-cell indicator. 

The original KXRCF shock detector is compared to the outlier-detection application 
in Figure [7| for /c = 2. The detected troubled cells at the final time using either the 
original or the outlier-detection approach are similar for k — 1. For A: = 2 and especially 
for A: = 3 fewer elements are detected by the outlier-detection scheme. However, more 
elements are detected in the top region of the domain. This is due to the fact that in 
this region neighboring jumps across the inflow edges of the element differ substantially 
from each other. 

The results using the minmod-based TVB indicator improve considerably when using 
outlier detection. In Figure [8(a/| the detected troubled cells at the final time are shown 
for the original minmod-based TVB indicator. Note that too many elements are detected: 
also continuous regions are selected. However, the outlier-detection technique applied to 
the DG coefficients only selects the correct discontinuity profile (Figure [8(b) [ ). It should 
be noticed that this approach detects discontinuities in the x- and ^/-direction, since DG 
coefficients are related to the first derivative in the x-direction, and to the 

first derivative in the ^/-direction. Fewer elements are detected in this case, and the 
approximation at time T = 0.2 is still accurate. 

5 Computational costs 

This section contains a discussion about the computational costs of the outlier-detection 
algorithm. First, we sort 2^“^ vectors of length 16 each. We use the ’Selection sort’ 
sorting algorithm, which finds the minimum value of the vector, swaps it with the value 
in the first position, and repeats these steps for the remainder of the list. The method 
is of order 0{N‘^) time complexity, but it is possible to use a more efficient sorting 
algorithm (for example of order 0{N)) [36]. Once the vectors are sorted, we compute 
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the quartiles and outer fences. Outliers are determined by comparing the smallest vector 
values with Qi — 3{Qs — Qi) and the biggest components with + 3{Qs — Qi). For 
the smallest values we start with testing whether dQ < Qi — 3{Qs — Qi). If dg 
outlier, then there are no other outliers, since di > do > Qi — 3{Qs — Qi). If d^ is an 
outlier, then we test df in the same way. Note that (by construction) at maximum dg, df 
and ^2 should be tested (as they are the only possible low outliers). Similarly we test 
df^ and (possibly) df^ and df^ against Qs + 3{Qs — Qi) (depending on the outcome). 
Finally, the detected outliers in the left half of the local regions are compared with the 
bounds of the left-neighboring region, and the outliers in the right half are compared 
with the right-neighboring region. 

It should be noticed that this novel method works well on a CPU. The local vectors 
can also be considered using parallel architectures. However, in that case the costs for 
communication will be higher, since local information should be distributed along the 
devices. On the other hand, it also typically results in fewer places where a limiter must 
be applied. 

In Table [1] the computational times are shown for the test problems of ® using 
either the original or the outlier-detection indication technique. Notice that the compu¬ 
tational times using outlier detection are slightly longer than the original times, except 
for the KXRCF indicator. In that case, the number of detected elements for the original 
algorithm is much larger than when outlier detection is applied, such that the moment 
limiter is applied more often. For the rest of the examples, the increase in computa¬ 
tional time is on average 2.9%, which is reasonable. It should be emphasized that the 
new method also reduces the number of tests by not having to find a problem-dependent 
parameter. 



Multiwavelets 

KXRCF 

Minmod 


Original 

Outlier 

Original 

Outlier 

Original 

Outlier 

Sod 

0.187 

0.208 

0.208 

0.212 

0.231 

0.256 

Lax 

0.263 

0.280 

0.299 

0.290 

0.329 

0.366 

blast wave 

10.539 

11.045 

13.505 

12.313 

14.776 

14.855 

Shu-Osher 

5.683 

5.845 

6.520 

6.512 

7.669 

7.973 


Table 1: Total computation time in seconds for the one-dimensional problems of 0 

The total computation times for the double Mach reflection problem are presented 
in Table [21 Note that the case k = 1 (minmod-based TVB indicator) is much faster than 
k — 2. Also here, the computation time increases, on average by 2.6%. Since no tests for 
parameter finding are needed, the new method will still provide the results much faster. 

6 Conclusion 

In this paper, we have introduced a new outlier-detection technique which can be applied 
to existing troubled-cell indication variables. In this way, problem-dependent parame¬ 
ters are no longer required. We showed the performance of this method for various 
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Multiwavelets 

KXRCF 

Minmod 

Original 

Outlier 

Original 

Outlier 

Original 

Outlier 

312 

316 

313 

324 

93 

97 


Table 2: Total computation time in minutes for the double Mach reflection problem 
{k — 2 for multiwavelet and KXRCF indicator, A: = 1 for minmod-based TVB indicator). 


test problems in one and two dimensions, using the modified multiwavelet troubled-cell 
indicator, the KXRCF shock detector, and the minmod-based TVB indicator. The re¬ 
sults were generally better than the original troubled-cell indicators using an optimized 
parameter: both the weak and the strong shock regions were detected, whereas smooth 
regions were not selected. Future work will be to include the local spatial information 
in the statistical approach, and to extend this to unstructured meshes. 
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(a) C = 0.05, a (b) Outlier, a 



(c) C = 0.05, 13 


(d) Outlier, jS 



(e) C = 0.05, 7 




(f) C = 0.05, total detected (g) Outlier, total detected 

Figure 6: Detected troubled cells at T = 0.2, double Mach reflection problem, modified 
multiwavelet troubled-cell indicator. Ax = A?/ = 1/128, k — 2. 
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(a) Original 


(b) Outlier 


Figure 7: Detected troubled cells at T = 0.2, double Mach reflection problem, KXRCF 
shock detector. Ax = Ay = 1/128, k = 2. 




(a) M = 100 


(b) Outlier 


Figure 8: Detected troubled cells at T = 0.2, double Mach reflection problem, minmod- 
based TVB indicator. Ax = Ay = 1/128, k = 1. 
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